function [gbest,gbestval,fitcount,fitness,gbestval_history]= PSO_func(fhd,Dimension,Particle_Number,Max_Gen,VRmin,VRmax,pos,bool_Plot,varargin)
%[gbest,gbestval,fitcount]= PSO_func('f8',3500,200000,30,30,-5.12,5.12)
rand('state',sum(100*clock));
me=Max_Gen;
ps=Particle_Number;
D=Dimension;
cc=[2 2];   %acceleration constants
iwt=0.9-(1:me).*(0.5./me);
% iwt=0.5.*ones(1,me);
if length(VRmin)==1
    VRmin=repmat(VRmin,1,D);
    VRmax=repmat(VRmax,1,D);
end
mv=0.5*(VRmax-VRmin);
VRmin=repmat(VRmin,ps,1);
VRmax=repmat(VRmax,ps,1);
Vmin=repmat(-mv,ps,1);
Vmax=-Vmin;
% pos=VRmin+(VRmax-VRmin).*rand(ps,D);

e=feval(fhd,pos',varargin{:});

fitcount=ps;
vel=Vmin+2.*Vmax.*rand(ps,D);%initialize the velocity of the particles
pbest=pos;
pbestval=e; %initialize the pbest and the pbest's fitness value
[gbestval,gbestid]=min(pbestval);
gworstval = max(e);
gbest=pbest(gbestid,:);%initialize the gbest and the gbest's fitness value
gbestrep=repmat(gbest,ps,1);
gbestval_history = zeros(fix(me/500),2);

if(bool_Plot)
    figure(1)
    haxes = plot( 0 , 0 );
    XArray = [1];
    YArray = [gbestval];
    title('PSO');
end
%     figure(1)
%     haxes = plot( 0 , 0 );
%     XArray = [];
%     YArray = [];
%     title('PSO');

for i=2:me        
    if(gbestval == gworstval)
        mass = repmat(1/ps,ps,1);
    else
        mass = (e'-gworstval)/(gbestval-gworstval);
    end
    mass = mass / sum(mass);
%     sort_mass = [mass (1:ps)'];
%     sort_mass = sortrows(sort_mass,1);
%     std_mass = std(mass);
%     if(mod(i,5)==0)
%         std_mass = std(mass);        
%         XArray = [ XArray i]; 
%         YArray = [ YArray sum(sort_mass(ps-fix(0.8*ps):ps,1))];
%         set( haxes , 'XData' , XArray , 'YData' , YArray );
%         drawnow 
%     end

    aa=cc(1).*rand(ps,D).*(pbest-pos)+cc(2).*rand(ps,D).*(gbestrep-pos);
%     vel=iwt(i).*vel+aa;
    vel=0.5.*vel+aa;
    vel=(vel>Vmax).*Vmax+(vel<=Vmax).*vel;
    vel=(vel<Vmin).*Vmin+(vel>=Vmin).*vel;
    pos=pos+vel;
    pos=((pos>=VRmin)&(pos<=VRmax)).*pos...
        +(pos<VRmin).*(VRmin+0.25.*(VRmax-VRmin).*rand(ps,D))+(pos>VRmax).*(VRmax-0.25.*(VRmax-VRmin).*rand(ps,D));
    e=feval(fhd,pos',varargin{:});
    fitcount=fitcount+ps;
    tmp=(pbestval<e);
    temp=repmat(tmp',1,D);
    pbest=temp.*pbest+(1-temp).*pos;
    pbestval=tmp.*pbestval+(1-tmp).*e;%update the pbest
    gworstval = max(e);
    [gbestval,tmp]=min(pbestval);
    gbest=pbest(tmp,:);
    gbestrep=repmat(gbest,ps,1);%update the gbest

    if(bool_Plot & mod(iter,10)== 0)
        XArray = [ XArray i]; 
        YArray = [ YArray gbestval];
        set( haxes , 'XData' , XArray , 'YData' , YArray );
        drawnow
    end
    if(mod(i,500) == 0)
        gbestval_history(fix(i/500),1) = gbestval;
        gbestval_history(fix(i/500),2) = i;
    end
end

fitness = e;
end


